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Summary 


The  paper  presents  a  numerical  simulation  for  mixed  mode  crack  initiation  based  on  the 
concepts  of  damage  mechanics.  A  model  with  two  scalar  damage  variables  is  introduced  for 
characterization  of  damage  in  a  material  element.  Then  a  tangent  modulus  tensor  is  derived  for 
damage-coupled  constitutive  equations.A  failure  criterion  is  developed  with  the  concept  of 
damage  accumulation  not  only  to  identify  the  location  of  damaged  element  with  the  crack 
initiation  angle  but  also  to  determine  the  critical  load  for  mixed  mode  fracture.  The  damage 
model  developed  is  incorporated  in  a  general-purpose  finite  element  program  ABAQUA  through 
its  UMAT  subroutine.  The  finite  element  program  is  then  used  to  perform  numerical  simulation 
for  pre-cracked  specimens  under  monotonic  tensile  loading.  The  thin  plates  are  made  of 
aluminum  alloy  and  particulate  composite  embedded  with  a  crack  of  inclined  angle  P=0°,  30°, 
45°  and  60°  for  mixed  mode  fracture  analysis.  The  predicted  crack  initiation  loads  and  the  angles 
of  crack  initiation  agree  well  with  the  test  results. 

Introduction 

Most  conventional  approaches  for  mixed  mode  fracture  prediction  are  based  on  the 
theory  of  fracture  mechanics.  The  fracture  parameters,  such  as  the  stress  intensity  factor  Kc  and 
the  strain  energy  release  rate  Gc  for  brittle  materials,  the  J-integral  and  COD  for  ductile 
materials,  have  been  widely  applied  to  conduct  fracture  analysis  of  engineering  components 
containing  a  macro-crack.  However,  an  important  mechanism  of  failure  in  the  materials  is 
attributed  to  the  presence  of  micro-cracks/voids.  These  micro-defects  result  in  changes  of 
mechanical  property  in  the  form  of  material  degradation  due  to  initiation,  growth  and 
coalescence  of  these  micro-defects.  Complete  avoidance  of  such  material  damage  is  not 
realistic,  especially  for  stress  analysis  at  a  crack  tip.  In  addition,  research  results  have  cast 
doubts  on  the  validity  of  the  J-integral  and  COD  as  intrinsic  material  properties  for  ductile 
fracture,  which  by  definition  should  be  independent  of  geometry  and  loading  history  [1-3]. 

The  material  damage  has  been  successfully  characterized  with  the  theory  of  damage 
mechanics  first  introduced  by  Kachanov  [4]  and  later  developed  by  many  researchers  [5,6], 
With  the  introduction  of  averaging  macro-variables  (damage  variables),  the  deterioration  of 
materials  as  a  result  of  the  nucleation  and  growth  of  distributed  material  micro-defects  can  be 
determined  quantitatively.  Chow  and  Wang  developed  an  anisotropic  model  which  has  been 
successfully  applied  to  analyze  mixed  mode  crack  initiation  and  propagation  of  aluminum  alloy 


2024-T3  [7-9].  The  model  is  based  on  a  second-order  damage  tensor.  Considerable  computing 
time  is  however  required  for  the  FEM  analysis  on  the  transformation  between  local  coordinate 
system  and  principal  damage  coordinate  system.  For  the  sake  of  computing  efficiency,  Chow 
and  Wei  have  recently  proposed  an  isotropic  damage  model  with  two  scalars  [10].  This  paper  is 
intended  to  present  an  investigation  on  the  application  of  the  proposed  damage  model  to 
characterize  ductile  behavior  of  two  materials,  including  damage-coupled  constitutive  equations, 
finite  element  formulation  and  ductile  fracture  of  mixed  mode  crack. 

Damage-Coupled  Constitutive  Equations 

The  gradual  deterioration  of  the  material  under  service  due  to  nucleation  and  growth  of 
micro-cracks  or  defects  can  be  characterized  with  an  internal  state  variable  known  as  damage 
variable  [5],  Chow  and  Wei  have  recently  developed  a  two-scalar  damage  model  to  evaluate 
damage  accumulation  under  both  monotonic  and  cyclic  loading  [10].  This  section  provides  a  brief 
description  of  the  model  required  in  the  following  sections  for  the  development  of  finite  element 
formulation. 

Following  the  damage  mechanics  theory,  the  effective  stress  is  defined  as 
a  =  M(D) :  cr  (1) 

where  a  is  the  true  stress  tensor,  D  is  the  damage  tensor  and  M(D)  is  the  damage  effect  tensor 


which  can  be  expressed  as 
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D  and  p  are  two  scalar  variables  to  characterize  damage  accumulation  in  materials.  The 
thermodynamic  conjugate  forces  of  the  damage  variables  D  and  p,  known  as  the  damage  energy 
release  rate,  can  be  derived  from  the  Helmholtz  free  energy  \\i  as 
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where  C  is  the  elastic  tensor  for  damaged  materials 
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The  tensor  Z  can  be  expressed , 

as 
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E  and  v  are  respectively  the  effective  Young's  modulus  and  effective  Poisson's  ratio  for  damaged 
material.  The  relationships  between  E  and  v  and  the  damage  variables  D  and  p  are  established  as 
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Eq  and  vq  are  the  values  of  Young's  modulus  and  Poisson's  ratio  for  intact  or  undamaged  material. 


The  elastic  law  of  damaged  material  can  be  expressed  in  the  effective  stress-effective  strain 


space  as 


a=C0  :se 


(7) 


where  Co  is  the  initial  elastic  tensor  for  undamaged  material.  The  yield  surface  is  postulated  with 
the  concept  of  the  effective  stress  as 


Fp(or,R)  =  crp- [R0  +  R(p)]  =  0 
where  ap  is  the  effective  equivalent  stress 
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ceq  is  the  Von-Mises  equivalent  stress,  Ho  is  the  plastic  characteristic  tensor  for  undamaged 
material,  R0  is  the  yield  stress,  p  is  the  effective  equivalent  plastic  strain,  and  R  is  the  strain 
hardening  threshold.  The  constitutive  equations  of  plasticity  for  damaged  materials  are 
accordingly  derived  in  the  effective  stress-effective  strain  space  as 


d  sp  -  X 


d£p 

p  da 


dp  =  Xp 


d(~R) 


(10) 


The  plastic  damage  surface  is  formulated  with  the  thermodynamic  conjugate  forces  of  the 
plastic  damage  variables  as 


Fj(Yd,B)  =  Yd-[Bo  +  B(w)]  =  0 


(11) 


where  is  the  equivalent  damage  energy  release  rate  postulated  as 
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(12) 


Bo  is  the  initial  plastic  damage  threshold,  B  is  the  plastic  damage  hardening,  w  is  the  overall  plastic- 
damage,  and  y  is  the  damage  evolution  coefficient.  The  plastic  damage  evolution  equations  are 


dD  =  -Xd 
dw  =  -Xd 


3Yd 

dFj_ 

SB 


=  4 


XjYd 

2Yd 


dM  =  -Xd^- 

dY„ 


XjTYj, 

2Yd 


where  Xd  is  the  Lagrange  multiplier. 


(13) 


Finite  Element  Formulation 

The  proposed  damage  model  is  discretized  and  coded  in  the  user  subroutine  UMAT  of  a 
finite  element  package  known  as  ABAQUS  (Version  5.8).  The  implementation  aims  at 
providing  a  tool  for  numerical  analysis  based  on  the  proposed  damage  model  and  validating  the 
model  by  comparing  its  predictions  with  experimental  measurements.  The  procedure  is  similar 
in  principle  to  the  conventional  FEM  analysis,  except  that  the  tangent  modulus  tensor  is  coupled 
with  the  damage  variables.  The  formulation  of  the  tangent  modulus  tensor  is  derived  as  follows. 

The  total  elastic  strain  shown  in  equation  (7)  can  be  alternatively  written  as 


Multiplying  the  above  equation  by 
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From  the  yield  surface  equation  (8),  we  have 
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Substituting  the  equations  (10)  and  (16)  into  (15),  we  obtain  the  incremental  plastic  strain  as 
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Thus  the  relationship  between  dcr  and  ds  is  obtained  with  equations  (10),  (14)  and  (17)  as 

da  =  Cp  :  ds  (18) 

where  Cep  is  the  instantaneous  tangent  modulus  tensor  expressed  as 
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From  the  definition  of  yield  surface  in  equation  (8), 
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Substituting  it  into  equation  (19),  the  tensor  CL0P  can  be  derived  with  equation  (9)  as 
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The  above  formulation  of  the  instantaneous  tangent  modulus  in  equation  (19)  or  (21)  is 
expressed  in  terms  of  the  incremental  effective  stress  and  strain  in  equation  (18).  However  a 
finite  element  program  such  as  ABAQUS  is  written  in  the  true  stress-true  strain  space. 
Therefore,  it  is  necessary  to  transform  equation  (18)  to  the  conventional  stress-strain  space  as: 


dcr  =  Cep  :de 


(22) 


where  Cep  is  the  effective  instantaneous  tangent  modulus  tensor  in  true  stress-true  strain  space 
which  can  be  derived  and  described  in  the  following  section. 

From  the  definition  of  effective  stress  in  equation  (1), 

da  =  dM:a  +  M:do-  (23) 

The  relationships  between  the  effective  and  true  elastic  strain  and  plastic  strain  are 
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Therefore,  the  effective  strain  is  derived  as 
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With  the  plastic  damage  evolution  laws  in  equation  (13),  dM  can  be  expressed  as 
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From  the  plastic  damage  surface  expressed  in  equation  (11) 
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Equation  (27)  becomes 
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Combining  equations  (18),  (22),  (23),  (25),  (26)  and  (29),  the  effective  instantaneous  tangent 
modulus  tensor  is  derived  instead  as 
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In  addition,  the  derivatives  in  equations  (29)  and  (31)  are  derived  as 
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Usually  the  matrix  form  of  Cepin  the  equation  (30)  is  a  6x6  asymmetric  matrix  for  which  most 
general  purpose  finite  element  programs  may  encounter  computational  convergency  difficulties. 
Consequently,  the  symmetric  form 


Cs  =  ^(Cep+  CpJ  ) 


(37) 


is  taken  as  a  tangent  modulus  matrix  when  the  finite  element  stiffness  matrix  is  computed 

K  =  JvB7C'B  (38) 


where  B  is  the  transformation  matrix. 

Crack  Initiation  Angles  for  A1 2024-T3  Plates 

A  thin  plate  of  aluminum  alloy  2024-T3  containing  an  isolated  crack  was  investigated 
first  as  shown  in  Fig.l.  The  length  dimension  L  is  86  mm,  the  thickness  of  the  plate  is  3.175  mm 
and  the  crack  length  is  1 5  mm.  0  is  the  inclined  angle  of  the  isolated  crack  and  p  is  the  angle 
measured  from  the  pre-crack  direction.  Different  values  of  0,  namely  0°,  30°,  45°  and  60°,  are 
chosen  for  numerical  simulation  with  ABAQUS  (Version  5.8).  Linear,  reduced-integration  solid 
elements  are  used  for  the  FE  analysis.  A  typical  finite  element  discretization  of  the  plate  with  a 


typical  inclined  crack  of  45°  is  depicted  in  Fig.  2.  Radial  elements  were  chosen  around  the  crack 
tip  for  the  convenience  of  determining  angular  distributions  of  overall  damage  at  the  crack  tip. 
The  mechanical  properties  of  AL024-T3,  which  were  determined  and  reported  in  reference  [10], 
are  summarized  as: 

E0=74300  MPa  v0=0.34  R0=330  MPa 

y=-0.4  Bo=0.936  MPa  wc=0.185 

The  overall  damage  distribution  in  the  plate  containing  a  typical  inclined  crack  of  45°  was 
calculated  as  shown  in  Fig.  3.  It  can  be  observed  from  the  figure  that  the  damage  accumulation 
is  confined  around  the  crack  tip  region.  Crack  initiation  is  postulated  to  occur  at  the  location 
when  the  overall  damage  reaches  its  critical  value.  The  direction  of  the  crack  extension  (3j  is 
determined  from  the  detailed  damage  distribution  at  the  crack  tip.  The  angular  distributions  of 
the  overall  damage  at  the  crack  tip  are  calculated  as  shown  in  Fig.  4  for  different  inclined  cracks 
of  0°,  30°,  45°  and  60°.  The  normalized  damage  at  a  constant  radial  distance  from  the  crack  tip 
was  plotted  against  the  angle  of  rotation  p.  The  location  of  unit  value  of  the  normalized  damage 
is  used  to  determine  the  direction  of  crack  extension  (i.e.  the  initiation  angle  Pp. 

Table  1  summarizes  the  predicted  crack  initiation  angles  for  the  mixed  mode  specimens. 
The  measured  and  predicted  results  based  on  an  anisotropic  damage  model  reported  by  Chow 
and  Wang  are  also  included  for  comparison  [8],  It  can  be  observed  from  the  table  that  both  of 
the  numerical  results  agree  well  with  the  measured  ones.  An  advantage  of  the  proposed  damage 
model  is  its  ease  of  computation  relative  to  the  anisotropic  damage  model  with  a  second-order 
tensor. 


Fracture  Analysis  for  Particulate  Composite  Plates 

The  boundary  between  the  particles  and  the  matrix  in  particulate  composites  constitutes 
the  source  of  micro-defects,  which  will  grow  and  coalescence  until  a  macro-crack  is  formed 
under  load.  To  ignore  these  micro-structural  changes  in  the  failure  prediction  by  the 
conventional  methods  is  not  therefore  considered  realistic.  The  material  deterioration  due  to 
damage  accumulation  should  be  taken  into  account  for  durability  analysis  with  the  theory  of 
damage  mechanics.  The  proposed  damage  model  is  accordingly  applied  to  predict  mixed  mode 
fracture  of  a  thin  plate  made  of  a  particulate  composite  as  shown  in  Fig.  1 .  The  dimensions  of  the 
plate  are  4  in.  x  4  in.  x  0.2  in.  The  length  of  pre-crack  is  1  in.  and  its  inclined  angle  is 
represented  by  0. 

In  order  to  determine  mechanical  properties  and  damage  parameters  for  the  particulate 
composite  material,  a  modified  dog-bone  specimen  is  chosen  for  the  measurement.  The 
specimen  is  incrementally  loaded  at  different  strains  until  final  rupture.  Upon  each  incremental 
maximum  strain,  the  specimen  is  unloaded  to  zero-stress,  forming  a  hysteretic  loop.  The 
effective  Young’s  modulus  is  determined  after  each  unloading.  The  change  of  effective  Young’s 
modulus  is  taken  as  material  degradation  to  evaluate  the  damage  variable  D  as  shown  in  Figs.  5 
and  6.  The  true  stress-true  strain  curve  of  the  material  is  depicted  in  Fig.  7.  The  value  of  DCi 
which  is  considered  an  intrinsic  material  property,  was  determined  to  be  0.3.  For  the  particulate 
composite,  the  value  of  Poisson’s  ratio  is  assumed  to  be  constant  during  loading  process. 


Accordingly,  both  the  damage  evolution  coefficient  y  and  the  damage  variable  p  are  considered 
insignificant  and  thus  assumed  to  be  zero.  Therefore,  the  damage  model  is  simplified  to  the 
conventional  isotropic  damage  model  with  one  scalar  damage  variable. 

The  FE  analysis  is  similar  to  the  case  example  of  AL2024-T3.  The  loading  is  applied 
through  the  displacement  boundary  condition.  The  element  is  postulated  to  be  fully  damaged 
and  therefore  failed  when  the  plastic  damage  D  at  its  integration  point  reaches  a  critical  value  Dc. 
The  load  applied  at  Dc  is  determined  as  crack  initiation  load  and  the  position  of  damaged 
integration  point  is  used  to  define  the  direction  of  crack  extension,  i.e.  the  crack  initiation  angle. 
Three  different  inclined  crack  angles,  0  =  0°,  30°  and  60°,  are  chosen  for  the  analysis.  The 
numerical  results  on  the  crack  initiation  angle  and  the  fracture  load  are  summarized  in  Tables  2 
and  3,  demonstrating  a  satisfactory  agreement  with  the  experimental  measurements. 

Summary 

Damage-coupled  finite  element  formulation  has  been  derived  for  the  proposed  damage 
model.  The  model  has  been  implemented  in  ABAQUS  (version  5.8)  through  its  UMAT 
subroutine  for  failure  analysis  of  mixed  mode  fracture.  A  failure  criterion  is  developed  to 
determine  the  crack  initiation  at  the  location  where  the  overall  damage  accumulation  reaches  a 
critical  value  of  the  material.  Therefore,  the  crack  initiation  angle  and  fracture  load  can  be 
determined  from  the  FEM  results.  Two  engineering  materials,  AL2024-T3  and  a  particulate 
composite,  are  chosen  for  the  investigation.  A  mixed  mode  fracture  analysis  was  carried  out  and 
the  predicted  results  are  compared  well  with  the  test  data. 
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Table  1  Crack  Initiation  Angle  (3,  for  A12024-T3 
(characteristic  crack  length:  0. 1  mm) 


inclined  angle  0 

test 

numerical  simulation 

proposed  model 

anisotropic  model 

0 

0 

0 

0 

30 

35.9 

37.5 

43 

45 

53.7 

52.5 

56 

60 

71.2 

67.5 

73 

Table  2  Crack  Initiation  Load  for  Particulate  Composite 
(characteristic  crack  length:  1 .4  mm) 


Pre-crack  angle 

0 

0 

30 

60 

Load 

(lb) 

prediction 

24.0 

27.2 

36.9 

test 

23.4 

27.0 

36.2 

Table  3  Crack  initiation  Angle  Pi  (°)  Particulate  Composite 
(characteristic  crack  length:  1 .4  mm) 


Pre-crack  angle  0 
(°) 

0 

30 

60 

Pi 

0 

prediction 

0 

28 

62 

test 

0 

33 

68 

SDV22  VALUE 


Fig.  3  Damage  distribution  contours  in  AL2024-T3  plate  for  0  =  45° 


Fig  4  Angular  distributions  of  damage  for  mixed  mode  AL2024-T3  specimen 
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Fig.  5  Normalized  Young’s  modulus  vs.  applied  strain  for  a  particulate  composite 


Fig.  6  Damage  evolution  curve  for  a  particulate  composite 


